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ABSTRACT 

We present results from a suite of Af-body simulations that follow the formation 
and accretion history of the terrestrial planets using a new parallel treecode that we 
have developed. We initially place 2000 equal size planetesimals between 0.5-4.0 
AU and the colli sional growth is followed until the completion of planetary accre- 
tion (> 100 Myr). A total of 64 simulations were carried out to explore sensitivity to 
the key parameters and initial conditions. All the important effect of gas in laminar 
disks are taken into account: the aerodynamic gas drag, the disk-planet interaction 
including Type I migration, and the global disk potential which causes inward mi- 
gration of secular resonances as the gas dissipates. We vary the initial total mass 
and spatial distribution of the planetesimals, the time scale of dissipation of nebular 
gas (which dissipates uniformly in space and exponentially in time), and orbits of 
Jupiter and Saturn. We end up with one to five planets in the terrestrial region. In 
order to maintain sufficient mass in this region in the presence of Type I migration, 
the time scale of gas dissipation needs to be 1-2 Myr. The final configurations and 
collisional histories strongly depend on the orbital eccentricity of Jupiter. If today's 
eccentricity of Jupiter is used, then most of bodies in the asteroidal region are swept 
up within the terrestrial region owing to the inward migration of the secular reso- 
nance, and giant impacts between protoplanets occur most commonly around 10 
Myr. If the orbital eccentricity of Jupiter is close to zero, as suggested in the Nice 
model, the effect of the secular resonance is negligible and a large amount of mass 
stays for a long period of time in the asteroidal region. With a circular orbit for 
Jupiter, giant impacts usually occur around 100 Myr, consistent with the accretion 
time scale indicated from isotope records. However, we inevitably have an Earth 
size planet at around 2 AU in this case. It is very difficult to obtain spatially concen- 
trated terrestrial planets together with very late giant impacts, as long as we include 
all the above effects of gas and assume initial disks similar to the minimum mass 
solar nebular. 

Key words: Accretion- Origin, Solar system- Planetary formation-Terrestrial 
planets 



1. Introduction 

Accretion models for the formation of the terrestrial planets need to reproduce the follow- 
ing properties of our solar system. (1) The spatial mass distribution of the terrestrial planets is 
radially concentrated. Namely, the Earth and Venus contain most of the mass inside of Jupiter's 
orbit, whereas the masses of Mars and Mercury are an order of magnitude smaller than the 
Earth's mass (M e ) and the total mass of asteroids is only 5 x 10 -4 M @ . The orbital separation 
between Earth and Venus is only 0.3 AU. (2) The orbits of the Earth and Venus are nearly cir- 
cular and coplanar even though orbital excitements are expected in the late stage of accretion 



where large protoplanets experience mutual collisions. This probably indicates that there were 
some damping mechanisms to reduce their orbital eccentricities and inclinations near the end of 
their formation processes. (3) Isotope records such as W/Hf and U/Pb for the Earth and Moon 
indicate that the Moon-forming giant impact occurred ~ 100 Myr after the beginning of the 
solar system (Touboul et al., 2007; Allegre et al., 2008). Some accretion models can reproduce 
two of the three properties but there is no accretion model that satisfies all three properties. 

/V-body simulations that follow only the gravitational interactions of protoplanets (or with 
the giant planets) usually produce too large orbital eccentricities of the final planets (Chambers, 
1998; Agnor et al., 1999; Raymond et al., 2004; Kokubo et al., 2006). To resolve this issue two 
damping mechanisms have been proposed: one is by remnant planetesimals (Chambers, 2001; 
O'Brien et al, 2006; Raymond et al, 2006, 2009; Morishima et al, 2008) and another is by 
remnant nebular gas (Agnor and Ward, 2002; Kominami et al., 2002, 2004; Nagasawa et al., 
2005; Ogihara et al, 2007; Thommes et al., 2008). 

In particular, the dynamical shake-up model of Nagasawa et al. (2005) and Thommes et 
al. (2008) is the most successful model so far in reproducing both the small orbital eccentric- 
ities and the radial mass concentration. Their model initially places a few tens of protoplanets 
embedded in a gas disk and takes into account the effect of the gas giant planets with orbital 
eccentricities comparable to or slightly larger than the present values (the present eccentrici- 
ties of both Jupiter and Saturn are ~ 0.05). The secular resonance with Jupiter migrates from 
~4.0 to 0.6 AU with dissipation of the gas disk and sweeps protoplanets toward the current 
location of the Earth. Mutual impacts between protoplanets occur even when a large amount 
of gas remains, and the remnant gas disk reduces the orbital eccentricities of formed Earth-like 
planets. 

If the formation of the Moon occurred very late, as indicated by recent isotope measure- 
ments (Touboul et al., 2007), significant damping due to remnant gas may be problematic for 
the following reasons. Significant depletion of gas already at ~ 10 Myr is suggested from di- 
rect measurements of line emissions of circumstellar disks (Pascucci et al. 2006) and from Ha 
emission profiles of central stars due to accretion shock (Muzerolle et al., 2000), although we 
cannot rule out the possibility that our protosolar nebular was exceptionally long-lived. Another 
potential problem for long-lived gas disks is depletion of the solid component due to Type I mi- 
gration, which is neglected in the model of Nagasawa et al. (2005) and Thommes et al. (2008). 
Once a protoplanet becomes as massive as Mars in a gas disk, the gravitational torque from spi- 
ral density waves launched by the protoplanet causes rapid inward migration of the protoplanet, 
referred to as Type I migration (Goldreich and Tremaine, 1980; Ward, 1986, 1997; Tanaka et 
al., 2002). McNeil et al. (2005) and Daisaka et al. (2006) conducted TY-body simulations taking 
Type I migration into account and found that the decay time scale for a gas disk needs to be ~ 1 
Myr in order to retain sufficient solid mass for the formation of Earth-size planets. 

Some remnant planetesimals are likely to have existed even at 100 Myr particularly in the 



asteroidal region. U Several authors have conducted A^-body simulations starting with proto- 
planets and remnant planetesimals (Chambers, 2001; Chambers and Cassen, 2002; O'Brien et 
al. 2006; Raymond et al. 2006, 2009). Their results show that final eccentricities of planets 
tend to be lower with increasing total mass of remnant planetesimals. Even with the same total 
mass of planetesimals, the damping is more effective when the size of an individual planetes- 
imal is smaller and planetesimals are placed uniformly over both the terrestrial and asteroidal 
regions rather than in the asteroidal region only (O'Brien et al., 2006; Raymond et al., 2006). 
An eccentricity of Jupiter comparable to or slightly larger than the present value is preferred in 
regard to the radial mass concentration (Chambers and Cassen, 2002; O'Brien et al., 2006; Ray- 
mond et al., 2009), although none of simulations have succeeded in reproducing the very high 
radial mass concentration of the current terrestrial planets. The spatial and size distributions 
of remnant planetesimals are arbitrarily assumed in these simulations. In order to clarify these 
distributions one needs to follow the growth of planetesimals and protoplanets in the runaway 
growth stage. In addition, the above simulations neglect the effects of gas. The protosolar neb- 
ular must have existed at least until the formation of Jupiter, and probably survived somewhat 
longer. Even a rapid dissipation of gas of the order of 0. 1 Myr in the presence of Jupiter strongly 
affects the orbital evolution of asteroidal bodies (Nagasawa et al., 2000). 

In the present study we perform iV-body simulations starting with a wide planetesimal 
disk and continue until the end of the accretion phase. Since we directly calculate the runaway 
growth of the protoplanets, we can obtain reasonable spatial and size distributions of remnant 
planetesimals in the late stage. We take into account the gas giant planets and all the possible 
effects of gas in laminar disks: aerodynamic gas drag, disk-planet interaction including Type 
I migration, and global disk potential. In a previous study, Morishima et al. (2008) also con- 
ducted Af-body simulations that followed the evolution of a protoplanetary disk starting with 
planetesimals only until the end of accretion. However, their simulations started from radially 
narrow annuli (widths of 0.3-0.5 AU) and neglected any external forces. In Sec. 2, we introduce 
our new Af-body code, review the effects of nebular gas, and outline our choice of simulation 
parameters. In Sec. 3, we show time evolution of two contrasting examples, present the results 
of all simulations and compare them with the solar system. In principle, none of our simula- 
tions succeed in satisfying all the three constraints listed in the first paragraph of this section. 
In Sec. 4, we discuss possible solutions for this problem. In Sec. 5, we give our conclusions. 



2. Method 

2.1. N-body code 

Any Af-body code consists of a gravity solver to calculate the mutual gravity force and an 
integrator to update positions and velocities for a next time step. Existing A^-body codes applied 



1 In the present paper, we loosely refer to the outside and inside of 2 AU as the asteroidal and terrestrial regions, 
respectively. 



for the runaway growth stage have good gravity solvers, such as tree methods (Richardson et al., 
2000; Stadel, 2001) or special hardware such as HARP or GRAPE (Kokubo and Ida, 1996). On 
the other hand, V-body codes applied for the late stage have good integrators, such as SyMBA 
(Duncan et al., 1998) and Mercury (Chambers, 1999), so that one can take a large time step 
size, while solving the mutual gravity by direct N 2 calculation. The new Af-body code we have 
developed has both a good gravity solver and a good integrator, allowing us to follow accretion 
stages of planetary systems with large N and over a large number of orbits. The code is still 
being upgraded and details of the complete version will be reported elsewhere. Here we briefly 
describe the version of the code used for simulations in this paper. 

The gravity solver we use is a parallel-tree method, PKDGRAV (Stadel, 2001). The orig- 
inal version of PKDGRAV uses a leap-frog integration scheme, where individual particles are 
updated in hierarchical time step blocks, each block with a factor of two smaller time-step 
according to the local dynamical time of its particles. Richardson et al. (2000) added some 
functions for planetary accretion to the code and conducted simulations for the runaway growth 
stage. In the tree method, all particles are divided into hierarchal tree-structured groups. Then, 
the gravity forces from distant particle groups are calculated using their multipole expansions. 
The zeroth-order multipole expansion corresponds to the gravity from the mass center of a par- 
ticle group. In order to keep good accuracy, we include multipole expansions up to the fourth 
order. The computer time of gravity calculation using this gravity solver is only proportional 
to N log N. We find a significant speed up in the gravity calculation with our tree method as 
compared with the direct N 2 calculation when N > 200. In addition, the tree method is suitable 
for parallel computing, and this results in further speed up when N > 1000. Since the number 
of particles decreases with growth of planets, we switch calculation methods from parallel-tree 
to serial-tree, and to serial-direct N 2 calculations. The error in the total energy due to the tree 
method appears as a random walk, and thus is proportional to the square root of the number of 
time steps. Therefore, with some short additional simulations, we can estimate the accumulated 
error after a long term evolution of a system. The accuracy of the method is controlled by the 
opening angle for tree cells and we set it to be 0.5 in the present paper. We confirm that the 
expected error due to the gravity calculation is sufficiently smaller than the error due to the 
integrator's, which rises most notably during close encounters. 

The integrator implemented in the code is a mixed- variable symplectic (MVS) integrator, 
SyMBA (Duncan et al., 1998). In MVS integrators such as SyMBA and Mercury, the Hamilto- 
nian is split into a part which describes the Keplerian motion around the central star and another 
part which describes the perturbation from other bodies. The former part corresponds to up- 
dates of positions and velocities along Kepler orbits. This can be done nearly analytically. The 
latter part corresponds to velocity changes due to the perturbation from other bodies, which we 
calculate using the tree method described above. These types of integrators are very accurate 
and allow us to take very large time steps for long term simulations. We use a time step size, 
At, of 6 days. This is comparable to those used in previous high-resolution V-body simulations 
using MVS integrators (O'Brien et al., 2006; Raymond et al., 2006). Practically, the time step 
size needs to be chosen so that a fraction of particles in close encounters defined below is small 
enough. 



In close encounter search, we first search for 10-20 nearest neighbors of each particle 
using the same algorithm in the tree SPH code, Gasoline (Wadsley et al., 2004). With the 
relative velocities and positions at the beginning and end of a time step, the closest distance 
of a pair during the step is calculated using a third order interpolation (Eq. (11) of Chambers, 
1999). In order to judge whether a particle is in a close encounter or not, we first define the 
critical distance for each particle. The critical distance is defined as the larger of (1) n\Ru or 
(2) H2V ep At, where n\ and n 2 are numerical factor (we use n\ = 4.5 and n 2 = 1.0), Ru is the 
Hill radius, and v ep = v^ sp Ve 2 + i 2 is the epicyclic velocity of the particle with v^ ep , e, and 
i being the Kepler velocity, eccentricity, and inclination, respectively. If the closest distance 
of a pair is smaller than the sum of the critical distances, we regard this pair to be in a close 
encounter. In the original SyMBA code (Duncan et al., 1998), only the first condition is used. 
However, the second condition is necessary if the relative velocity is much larger than the mutual 
escape velocity; such encounters usually occur near resonances of the giant planets or between 
small planetesimals stirred by large protoplanets. In the Mercury code (Chambers, 1999), the 
second condition uses the largest v kep in all particles instead of each body's v ep which we apply. 
Although this criterion of the Mercury code gives better energy conservation (for a same n 2 ), 
it also results in too many pairs being identified as involved in close encounters during the 
runaway growth stage, where the number density of particles is large. If the number density 
is not large, as is the case in the late stage, use of the Mercury's criterion is recommended 
(although we do not use it even in the late stage in this paper). 

All the above computations are performed in parallel (if N > 1000), while particles in 
close encounters are collected on a single cpu during the domain decomposition; the procedure 
which distributes the workload, and hence the remaining particles, amongst processors. Then, 
multistep time stepping is done until all particles proceed to the end of the original time step, 
following the procedure described in Duncan et al. (1998). The closest mutual distance in a 
substep is always estimated by interpolation using the information of both the beginning and 
the end of the step. This guarantees time symmetry of the integrator which is essential to avoid 
secular errors. The computational time for the substep procedure is small enough as it is simply 
proportional to the number of particles in close encounters. We set the ratio of the cut-off radius 
to that of the next order to be 2.08 and the number of next order substeps in the current substep 
to be 3 after Duncan et al. (1998). 

We assume that two bodies merge when they partially overlap. We assume the density of 
particles to be 2 g cm 3 . Since we do not apply the artificial enhancement of radii (Kokubo 
and Ida, 2002, Morishima et al., 2008), the accretion time scale we will show is realistic. A 
particle is assumed to collide with the Sun when either the heliocentric distance is less than 
0.1 AU or the semimajor axis is less than 0.25 AU. A particle is assumed to escape from the 
system when the heliocentric distance is larger than 100 AU. The energy error after completion 
of accretion is typically AE/E = 10" 5 -10 3 . Note that energy changes due to merging or gas 
forces can be explicitly calculated and are subtracted in the above error. The largest portion of 
the error usually accumulates after the runaway growth stage where a large number of excited 
small particles experience mutual encounters. Therefore, the energy error primarily depends on 
how long small particles survive. 



2.2. Effects of gas 

In this section, we review the forces due to the nebular gas. The acceleration due to the gas 
is incorporated in the integrator by placing the half kicks before and after the Kepler drift oper- 
ation. This symmetric placement once again avoids secular errors coming from the inclusion of 
gas forces. 



2.2.1. Aerodynamic drag 

The aerodynamic drag force per unit mass (or acceleration) is given by (Adachi et al., 

1976) 

1 2 

/drag = ~Y~ C D7"- p p gas |i; rel |l> re i, (1) 

where m and r p are the mass and the radius of a particle, cd is the numerical coefficient (assumed 
to be 2 in the present paper), p gas is the gas density, and » re i = v - w gas is the velocity of a 
particle relative to the gas velocity. In the cylindrical coordinate (r, 6, z), the gas velocity in the 
9 component is obtained by the force balance in the r component as 

v 2 gas (r,z) c 2 / din Pgas dlnT\ GM Q r , N 

+ ^n — +71 ^F r -fgiobAr,z), (2) 



r r \ d In r d In r ) (r 2 + z 2 ) 3/2 

where c is the isothermal sound velocity, T is the gas temperature, G is the gravitational con- 
stant, and M Q is the mass of the central star (assumed to be the solar mass). The left hand side of 
Eq. © represents the centrifugal force while the first, second, and third terms of the right hand 
side represent the pressure gradient, the gravity of the central star, and the self- gravity of the gas 
disk (Eq. ([TBI)), respectively. Assuming the gas is isothermal in the z direction and neglecting 
the self gravity of the gas, the gas density is approximately given by 

ft " (r ' d = S exp (-|^ (3) 

where £ gas is the surface density of the gas and h = c/f2 kep is the scale height, with the Keplerian 
frequency Qkep- If we assume 2 gas oc r~ p , T oc r~ q , we obtain, 



dlnp gas tq 3\i, ( z 



amr =-" + [i-2J[ l -[- h ))- (4) 

We adopt the temperature profile of an optically thin disk, T = 280 (r/1 AlTr 1/2 K (q = 1/2), 
after Hayashi et al. (1981). We will explain how S gas and p are given in Sec. 2.3. 

The growth rate of a protoplanet is proportional to e~ 2 , where e is the orbital eccentricity 
of planetesimals, and the balance between viscous stirring by protoplanets and damping due to 
gas drag gives e oc m 1/18 (Eq. (25) of Kokubo and Ida, 2000). The initial planetesimal mass 
m in our simulations is ~ 10 25 g. On the other hand, the size of planetesimals m G i formed by 



gravitational instability is estimated to be several order of magnitude smaller (Goldreich and 
Ward, 1973). Therefore, the growth time scale of protoplanets will be significantly overesti- 
mated in our model. In order to reproduce the growth rate of protoplanets surrounded by small 
planetesimals (~ m G i), we artificially enhance the gas drag force on planetesimals by a factor 
of (m/m') 2/3 in /V-body simulations if the particle mass m is smaller than ni\ (> m ). Here the 
mass m! is given as 

log m' - log m Gl log m - log m 

= . (for m < mi) (5) 

log m,\ — log niQi log nil — log mo 

With this scheme, a bunch of small planetesimals of mass m' are represented by a single super 
planetesimal of mass m (e.g., McNeil et al., 2005). The mass nil should be as small as possible 
in order to avoid artificial clean up of large planetesimals by protoplanets. We use ni\ = 0.01 M 
as a compromise, and m GI = 10 19 g. The migration and damping time scales with the artificial 
enhancement of the gas drag are shown in Fig. 1 for the case of e = 0.05 and m = 0.0025 M m . [Fig. 1] 



2.2.2. Tidal interaction between a protoplanet and a nebular disk 

We use the formulation of Tanaka et al. (2002) for Type I migration rate. For the damping 
rates of orbital eccentricity e and inclination i, we use the formulation given in Tanaka and Ward 
(2004). We make corrections in these formulations for the case of large e and i after Papaloizou 
and Larwood (2000) and Cresswell et al. (2007). 

Tanaka et al. (2002) derived the torque due to tidal interaction between a protoplanet and 
a three-dimensional baroclinic disk, for the case of e = i = 0. The decay time scale of the 
semimajor axis a is given as 

a . 1 ^ke P ,„,,,.,., N -i M © M ° / c \ 2 i 



rtidi = -t = -t-^ = (2.7 + \A P r-^—^ — Q k ;, (6) 

where / t idi,6> is the acceleration in the azimuthal direction. 

Tanaka and Ward (2004) extended their work to the case of nonzero e and / (<sc hi a) and 
obtained the rate for damping of e and i due to the tidal interaction. The r, 6, and z components 
of the damping force per unit mass are given as (a typo was corrected; see Ogihara et al., 2007) 

fm, r = [0.104[^-4 p (r)]+0.176iv], (7) 

^"wave 

fufi = —[-1.736^ -i4 5p (r)]+0.325i V ], (8) 

'wave 

/ticu = — [- 1.088^ + 0.87 lzQ k ep(r)], (9) 

~wave 

where characteristic time of the orbital evolution, r wave , is given by 

1 /V „2\-l / „ \4 
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In Eqs. ©-©, the velocity components of a particle are given as v = (v r , vg, v z ) and v' ke (r) = 
^L _ r fgiobA r ' O)) 1 '' 2 is the Keplerian velocity of solid bodies including the gravity of the gas 
disk. Note that all the components of the damping force become zero when e = i = 0. 

The formulation of Tanaka and Ward (2004) is limited to small e and i («c h/a), or that 
the planet motion relative to the gas remains subsonic. From a two dimensional linear theory, 
Papaloizou and Larwood (2000) found that the time scales r t idi and r wave (Eqs. © and (PTOl) ) 
increase with increasing e by factors of g t idi and g wave , respectively: 

_ l + Cx/1-3) 5 

gtm -i-( X /i.ir ( } 

i , 

#wave = l + ~X . (12) 

Here the original x used in Papaloizou and Larwood (2000) is x PL = ae/h. On the other hand, 
we use a modified value as 

X= 2h ' (13) 

in order to reproduce similar results of simulations of Cresswell et al. (2007) (see the next 
paragraph). The evolution of the orbital eccentricity with the corrections for large e is given as 

m\ = _^^. (14) 

c \"'/tjd ' wave y wave 

This equation corresponds to Eq. (45) of Tanaka and Ward (2004), but the correction factor is 
multiplied by T wave . From Eqs. ©, (fTT|) . and fp~4|) . the evolution of the semimajor axis due to 
the tidal interaction is given by, 

1 / da\ 1 1 / de\ / 2e 2 , 



a \dt Aid T tidi ^tidi e\dt ) tid \ Vi _ t 

where the first term in the right hand side comes from exchange of orbital angular momenta 
while the second term represents the migration rate due to energy dissipation. 

Cresswell et al. (2007) conducted hydrodynamic simulations for migration of a 20 Earth 
mass planet with nonzero e and i embedded in a gas disk. They found good agreement with 
Tanaka et al. (2002) and Tanaka and Ward (2004) for small e and / (< 0.1), while the damping 
time scale e/e is proportional to e 3 for large e as predicted by Papaloizou and Larwood (2000) 
(Eq. CG]))- They also showed the same dependence of the inclination decay rate on i, so we 
first simply replace e by (e 2 + i 2 ) 112 in Eq. (fT3l) . With the original jcpl, Eq. (fT5l) gives outward 
migration for ae/h > 1.1. Cresswell et al. (2007) also found that the sign of the torque changes 
to positive (g^a < 0) when ae/h > 1.6, but did not find any indication of outward migration. 
This means that the absolute value of the second term in the right hand side of Eq. (TT5T) is 
always larger than the first term for large e. It is unclear whether outward migration really 
occurs or not and we need to wait for future studies. In the present paper, we use the modified x 
(Eq. (TT3T) ) to reproduce similar dependences of e and a on e and i to those obtained in Cresswell 
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et al. (2007). The migration and damping time scales (e/e and a/a) with various e are shown 
in Fig. 1. Although there are some differences between the time scales shown in Fig. 1 and 
those from Cresswell et al. (2007), these do not result in significant differences in our /V-body 
simulations unless the time scales become too large or negative. Practically, large protoplanets 
obtain large e when they are trapped in secular resonances. If outward migration really occurs, 
protoplanets can be removed from secular resonances which migrate inward. 



2.2.3. Global nebular force and secular resonances 

We calculate the gravity of the gas disk after Nagasawa et al. (2000). The r and z compo- 
nents of the nebular gravitational force are given as 



/glob,r(7\ Z) 



- 2G ffv 



Pgas(r',Z') 



X 



V(r + r'Y + (z - Z') 2 
r 2 - r' 2 - (z - z') 2 



fglob,z( r ' z ) 



-E(k) + K(k) 

(r - r'Y - (z - zY 

p g , s (r',z')r'(z-z') 



dr'dz 



-4G r r 

U -co \J r' 



V(r + r'Y + (z - z') 2 



E(k) 



(r - r'Y + (Z- z') 



'\2 



dr'dz, (16) 



where K(k) and E(k) are elliptic integrals of the first and second kind, and k is given by 



k = 



Arr' 



(r + r'Y + (z- z'Y 



(17) 



These components of the force are initially tabulated in a (r, z) grid, and interpolated in sim- 
ulations. A secular resonance occurs when precession rates of two orbiting bodies coincide 
(Heppenheimer, 1980; Ward, 1981; Nagasawa et al., 2000). We analytically calculate preces- 
sion rates of planetesimals and the gas giant planets for e, i <sc 1 , following Appendix A of 
Nagasawa et al. (2000). Then, we obtain the locations of the secular resonances with Jupiter 
(v 5 ) and Saturn (v 6 ) (e.g., Fig. 2). 



2.3. Input parameters 



We adopt four different initial spatial distributions of planetesimals, four different dissipa- 
tion time scales of the nebular gas, and four different types of orbits of Jovian planets. In total, 
therefore, we made 64 runs with different combinations of parameters. Each run takes one to 
three cpu months until completion of accretion. 
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2.3.1. Initial spatial distribution of planetesimals 

We initially place 2000 equal-mass planetesimals between 0.5 and 4.0 AU. The initial 
surface density distribution of planetesimals is given by 

Ssoiid = £ solid , (y-^-j) P , (for 0.5 AU < r < 4.0 AU) (18) 

where £ S oiid,o is £ so iid at 1 AU. We adopt p = 1 and 2 in simulations. The initial total solid mass 
is given as 

£.0 AU 
2nr^ oM dr. (19) 

.. AU 

We adopt M T = 5 and 10 M e . In the case of M T = 5 M 9 , S so iid,o = 6.1 and 10.2 gem 2 for p = 1 
and 2, respectively. For comparison, the minimum mass solar nebula (MMSN) by Hayashi et 
al. (1981) has S solid , = 7.1 gem" 2 , p = 1.5, and M T = 4.3 M 9 . 



2.3.2. Dissipation time scale of nebular gas 

In the present paper, we limit ourselves to the case where the surface density of gas dissi- 
pates exponentially in time and uniformly in space: 



2. 



sir, t) = 2 gas , (y-^j) " exp (-^—\ (for 0.1 AU < r < 36.0 AU) (20) 

where 2 gas is 2 gas at 1AU and t = and Td eC ay is the time scale for gas decay. We adopt four 
different decay times: Td eC ay = 1, 2, 3, and 5 Myr. We assume that £ gaSi o = 2000 and 3366 gem" 2 
for p = 1 and 2, respectively. The gas-to-solid ratio is constant between 0.5 AU and 4.0 AU 
in the beginning of simulations and is 329 and 165 for M T = 5 and 10 M e , respectively. For 
comparison, MMSN has £ gaSj o = 1700 gem" 2 (Hayashi et al., 1981). We do not consider a gap 
around Jupiter's orbit in the present paper. The formation of a gap does not significantly alter 
the time of the sweeping secular resonance's passage in the gaseous region (Nagasawa et al., 
2005). However, orbital evolutions of asteroidal bodies in a gas-free gap will be different, in 
particular, if a gap is very wide. 



2.3.3. Orbits of Jovian planets 

We introduce Jupiter and Saturn after Tdecay from the beginning of a simulation, assuming 
that they form via core accretion (Bodenheimer and Polack, 1986; Dtoma et al., 2000; Lissauer 
et al., 2009). Rapid formation of Jovian planets via disk instability is unlikely for our solar 
system because efficient cooling for clump formation is only possible at r > 100 AU (Boley, 
2009). 
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We define aj and ej to be the semimajor axis and the orbital eccentricity of Jupiter, and 
a s and e s to be the corresponding orbital elements for Saturn. We adopt four types of orbits 
for Jupiter and Saturn at their formation and use the same notation introduced in O'Brien et al. 
(2006) and Raymond et al. (2009), although we have slightly different orbits: 

• EJS ("Eccentric Jupiter and Saturn"). Jupiter and Saturn are placed on their current orbits: 
a : = 5.20 AU, e : = 0.048, a s = 9.55 AU, and e s = 0.056. 

• CJS ("Circular Jupiter and Saturn"). These are the initial conditions used in the Nice 
model (Tsiganis et al., 2005). Jupiter and Saturn are radially more confined with circular 
orbits: aj = 5.45 AU, ej = 0.0, as = 8.18 AU, and e s = 0.0. The mutual inclination is 0.5 
degrees. 

• EEJS ("Extra Eccentric Jupiter and Saturn"). The same as EJS except a higher eccentric- 
ity of Jupiter: e } = 0.1. 

• CJSECC. Jupiter and Saturn have the CJS semimajor axes with higher eccentricities: 
e : = 0.05 and e s = 0.05. 

The EJS and EEJS simulations assume that Jupiter and Saturn formed at their current 
locations. On the other hand, the Nice model suggests that Jupiter migrated inward and Sat- 
urn migrated outward from their original CJS orbits via interactions with Kuiper Belt objects 
leading to a sudden change to the EJS orbits when Jupiter and Saturn crossed their 1:2 orbital 
resonance (Tsiganis et al., 2005). Asteroidal and cometary bodies strongly perturbed by this 
event are considered to have caused the Late Heavy Bombardment at ~ 3.85 billion years ago; 
well after formation of the terrestrial planets (Gomes et al., 2005; Strom et al., 2005). Initially 
higher ej and e s in the EEJS and CJSECC simulations than those in the EJS and CJS simula- 
tions, respectively, are considered because <?j and e s decrease via interaction with planetesimals 
(Sec. 3). 

We only take into account the gravity of other bodies and the global nebular force on 
motion of Jupiter and Saturn. We neglect any further interaction between a gas disk and the 
giant planets. Details of orbital evolution of a giant planet, which opens a gap in a gas disk, have 
a complicated dependence on disk properties, such as shape of a gap and disk mass, for which 
explicit and reliable formulae are not yet available. Nevertheless, some basic trends, which 
should be included in future studies, started to be clarified as follows. Analytic (Goldreich and 
Sari, 2003; Ogilvie and Lubow, 2003) and numerical (D'Angelo et al., 2006) studies indicate 
that eccentricity growth for a Jupiter-mass planet occurs together with the saturation of the 
corotation resonance in a low-viscosity disk. The migration time scale of a Jovian planet is 
much longer than the viscous diffusion time scale and even outward migration occurs, either 
when the gap is clean and the orbit is eccentric (D'Angelo et al., 2006) or when the gap is 
not deep enough (Crida and Morbidelli, 2007). A long migration time scale of Jupiter due to 
disk-sculpturing by Saturn is also pointed out by Morbidelli and Crida (2007). 
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3. Results 
3.1. Time evolution 

3.1.1. Case with EJS orbits 

Figures 2-4 show an example of the time evolution for the case with the present orbits of 
Jupiter and Saturn: snapshots on the plane of the semimajor axis versus the orbital eccentricity 
(Fig. 2), evolution of the total masses and the mean eccentricities of small and big particles 
(Fig. 3), and evolution of the masses and orbital excursions of surviving planets (Fig. 4). The gas [Figs. 
decay time is t& ecay = 1 Myr, the initial total mass is 5 M 9 , the power-law index for the surface 4] 
density is p = 2. In the present paper, we distinguish between a big particle, or (proto)planet, 
defined as a body with m > 2 x 10 26 g (~ 0.03 M e ), and smaller bodies or planetesimals. This 
border mass is slightly smaller than Mercury's mass (3.3 x 10 26 g), and one order of magnitude 
smaller than the isolation mass at 1AU (Kokubo and Ida, 2002). 

We follow growth of protoplanets without Jupiter and Saturn, until ? decay from the begin- 
ning of the simulation. We first observe the mass loss (~ 0.2 M & at 10 5 yr; see Fig. 3) by the 
migration of the smallest planetesimals due to gas drag at the inner most region. Protoplanets 
grow rapidly in the inner regions while many smaller bodies remain unaccreted in the outer 
regions. Protoplanets first tend to collide with smaller planetesimals because their orbital ec- 
centricities are smaller. At 1 Myr, the smallest particles significantly deplete at r < 2 AU, while 
some larger planetesimals still survive there (Fig. 2). Once a planet reaches a certain size, it 
starts rapid inward migration (Type I migration). Since the growth time is shorter for smaller r, 
significant mass depletion first occurs in the inner disk. With increasing time, the mass deple- 
tion propagates outward (see Ida and Lin (2008) for more details). The mass loss due to Type 
I migration is easily recognized in the upper panel of Fig. 3, as we see stepwise depletion in 
the total mass. More than half of the total mass is lost by this mechanism until 3 Myr. The 
half of the total mass corresponds to the initial total mass inside of ~ 1 .5 AU. The three planets, 
which survive till the end, are only lunar-size at 1 Myr and located outside of 1.2 AU (Fig. 4). 
There are other larger bodies located inside of 1.2 AU at 1 Myr, but they will not survive Type 
I migration. 

At 1 Myr, we introduce Jupiter and Saturn with their present orbital elements. Their for- 
mation causes a sudden increase of orbital eccentricities of nearby particles, particularly those 
in mean and secular resonances, while such a jump does not appear for inner planets (see the 
lower panel of Fig. 3). The locations of the v 5 and v^ resonances at t = 1 Myr are 3.10 and 3.34 
AU, respectively. These resonances, in particular v 5 , effectively enhance orbital eccentricities 
of planetesimals and protoplanets, and these bodies migrate inward as the gas drag and the tidal 
interaction damp their eccentricities. At 3-5 Myr, the mean eccentricities of big and small bod- 
ies are comparable (Fig. 3) as many protoplanets are also involved in the v 5 resonance. Near 
the resonance, the surface density of solid material is highly enhanced and relative velocities 
of solid bodies are small due to orbital alignments (Nagasawa et al., 2005). Consequently, pro- 
toplanets grow very quickly (Fig. 4). Collisions between protoplanets also commonly occur at 
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this stage. 

When the gas density becomes so low that inward migration of planetesimals and proto- 
planets due to the gas drag or the tidal interaction is slower than the migration rate of the secular 
resonance, these bodies are no longer swept by the resonance (Nagasawa et al., 2005). After the 
passage of the v 5 resonance, orbital eccentricities of planets are reduced first by tidal interaction 
with a remnant gas disk. The damping due to gas is effective as long as the damping time is 
shorter than the gas decay time (r damp < Td eca y)- F° r an Earth-size planet at 1 AU, this condition 
is satisfied until ~ 7-8 Tdecay Indeed, we see a rapid decrease of the mean eccentricity of big 
bodies at 5-7 Myr in Fig. 3. Later, the damping due to dynamical friction of planetesimals fur- 
ther reduces the eccentricities of planets. Its effect is determined by the total mass and the mass 
distribution of planetesimals (Morishima et al., 2008). Small bodies contain ~ 20% of the total 
mass at 10 Myr and this amount is likely to be sufficient to achieve energy equipartition between 
big and small bodies as long as the planet-planet interaction is negligible. The typical size of 
remnant planetesimals is the lunar size or slightly smaller. There are only a few remnant plan- 
etesimals with the initial size because the smallest planetesimals had small orbital eccentricities 
due to strong gas drag and efficiently merged with larger bodies. The equilibrium eccentricity 
of Earth-size bodies due to dynamical friction of lunar-size planetesimals is estimated to be ~ 
0.03 (Morishima et al., 2008), and we obtain a quite consistent value after most of the remnant 
planetesimals merge with planets. 

Through the above processes, three planets form near 1 AU at the end of this simulation. 
The largest planet has 1.4 M B while both of other two have 0.3 M @ . The locations of these 
three planets resemble those of the current terrestrial planets (except we do not have a Mercury 
analog), although the masses are somewhat different. Orbital eccentricities are as small as 
those for the current planets. The most serious issue for this simulation is that the last of the 
potentially Moon-forming impacts occurs too early (at 13.5 Myr; see Sec. 3.2 for the definition) 
when compared with that suggested from isotope records (Touboul et al., 2007). 



3.1.2. Case with CJS orbits 

Figures 5-7 show the time evolution for the case with circular orbits of Jupiter and Saturn, 
which are given the initial configuration of the Nice model. Except for the orbits of the gas [Figs. 5- 
giant planets, all other parameters are the same as those used in the EJS simulation in Sec. 3.1.1 7] 
(Figs. 2-4). If the orbits of the gas giant planets are nearly circular, the effects of the secular 
resonances are much weaker than the case of the EJS orbits. Consequently, a large fraction 
of particles stay in the asteroidal region for a long period of time. Nevertheless, the orbital 
eccentricities of the gas giant planets can not be completely zero when there are two (or more) 
gas giant planets. Therefore, resonances gradually sculpture orbits of particles in the asteroidal 
region. Some enhancements at the 3:1 mean motion resonance (2.6 AU) and the v 6 secular 
resonance (3.2 AU) are seen at 3 Myr (Fig. 5). These enhancements gradually propagate to 
non-resonant locations via mutual gravitational interactions (Raymond et al., 2006). The mass 
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transfer from the asteroidal region to the terrestrial region slowly occurs as excited asteroidal 
bodies encounter and collide with inner planets. 

At 3 Myr, there are about 15 Mars-size protoplanets in the terrestrial region. The fraction of 
small particles is very small in the terrestrial region whereas small particles are still dominant in 
mass in the asteroidal region. This situation is similar to some of the initial conditions adopted 
in previous /V-body simulations (e.g., Raymond et al., 2006). In our simulation, however, the 
mean orbital separation between protoplanets is ~ 20 Hill radius at 3 Myr while previous late 
stage simulations usually assume the orbital separation to be ~ 10 Hill radius, based on the 
oligarchic growth model of Kokubo and Ida (1998). The larger separation in our simulation 
is due to Type I migration, which makes inner protoplanets start to migrate inward earlier than 
outer planets (McNeil et al., 2005). With the large orbital separation and the orbital stabilization 
due to the nebular gas, mutual collisions between protoplanets are rare until ~ 10 Myr. 

After 10 Myr, mutual collisions between protoplanets start to occur, and the orbital separa- 
tion increases as the number of protoplanets decreases. Strictly speaking, the orbital separation 
normalized by the mutual Hill radius, and averaged over pairs of planets, increases nearly log- 
arithmically with time (from ~ 20 at 3 Myr to ~ 60 at 800 Myr), and this is probably because 
the stable time scale of a planetary system exponentially increases with the normalized separa- 
tion (Chambers et al., 1996; Ito and Tanikawa, 1999; Yoshinaga et al., 1999). Through giant 
impacts, an Earth-size planet forms slightly inside of 2 AU. This planet highly enhances orbital 
eccentricities of asteroidal bodies and their number gradually decreases by either hitting this 
planet or escaping from the system after encounters with Jupiter. 

There are still a large number of remnant planetesimals at 100 Myr, but not at 300 Myr 
(Figs. 5 and 6). Therefore, damping due to planetesimals is expected until giant impacts com- 
plete at around 100 Myr. In this simulation, the last potentially Moon-forming impact occurs at 
139.5 Myr on the outermost planet (Fig. 7). Nevertheless, the mean orbital eccentricity in the 
end of the simulation (~ 800 Myr) is not small, as the two middle planets are still strongly inter- 
acting (Fig. 7). These two planets are stirred by the innermost and outermost planets which have 
larger masses than the middle planets. This simulation succeeded in reproducing the very late 
giant impact, but failed to reproduce the nearly circular orbits and the radial mass concentration 
of the current system. 



3.2. Comparison with constraints 

In this section, we examine the properties of planetary systems formed in all runs and 
compare them with those for our own system. We consider the following three properties of 
systems (see also Raymond et al., 2009). 

1. The spatial mass distribution. The radial mass concentration statistic, S c , is defined as 
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(Chambers, 2001) 

S c = max|- =-^ |, (21) 

\Lmj\og l0 (a/aj)) 

where m 7 and cij are the mass and the semimajor axis of each planet, and the subscript j is 
given in the order of a,-. The function in the parenthesis is calculated for a throughout the 
terrestrial planet region, with S c being the maximum value. This quantity increases with 
radial mass concentration. We use two additional parameters as constraints. The first one 
is the minimum orbital separation between neighboring planets b m [ n normalized by their 
mutual Hill radius r H j = 0.5(a ; - + a/+i)[(»i/ + m i+1 )/(3M )] 1/3 : 

b min = mm . (22) 

\ r H,y / 

The second being the mass weighted mean semimajor axis of planets 

Z am 



M- 



(23) 



T.final 



where M T> fi na i = 2 m j is the final total mass of planets. For the terrestrial planets, S c = 
89.9, b mm = 26.3 (the normalized orbital separation between the Earth and Venus), and 
a m =0.90 AU. 

2. The deviation from circular and coplanar orbits. The orbital excitation of terrestrial plan- 
ets is characterized by the angular momentum deficit, S d (Laskar, 1997): 



S d = 



Y,mj^a~j]A- cos0' 7 ) ^1 - efj 



Zmj^a] 



2\\ 

(24) 



where e 7 and ij are the orbital eccentricity and inclination of each planet with respect to 
the fiducial plane of an initial planetesimal disk. We also measure the mass weighted 
mean eccentricity: 

g e J m J n ,, 

£m = M • (25) 

M T,final 

We use the orbital elements averaged over a few Myr. For the terrestrial planets, S d = 
0.0018 and e m = 0.038 (Quinn et al., 1991). 

3. The timing of the Moon-forming impact. Based on simulations for giant impacts and 
accretion of moons (Cameron and Benz, 1991; Idaetal., 1997; Canupetal., 2001; Canup, 
2004; 2008), we roughly regard an impact as potentially Moon-forming, if a) the total 
mass of the impactor and the target is > 0.5 M @ , b) the impactor's mass (< the target's 
mass) is > 0.05 M m , and c) the impact angular momentum is > 0.5 L EM , where L EM is the 
current total angular momentum of the Earth-Moon system (i.e., the Earth's spin angular 
momentum and the Moon's orbital angular momentum). We define the time of the last 
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potentially Moon-forming impact in a run to be t{ mp . The time of the actual Moon-formin 
impact is estimated as t imp = 50-150 Myr (Touboul et al., 2007; Allegre et al., 2008). 
The small content of iron in the Moon is explained when the impact velocity is as slow 
as the mutual escape velocity (Canup, 2004, 2008); this constraint will be also discussed 
(Figs. 13 and 16). 

Raymond et al. (2009) considered two additional constraints. The first one is the absence 
of planetary bodies (> 0.05 M B ) in the asteroidal region (> 2 AU). They showed that bodies 
which are too massive destroy the observed radial variation of the asteroidal taxonomic types 
(Gradie and Tedesco, 1982). If massive bodies exist in the asteroidal region, this system usually 
has too large a m or too small S c , because we do not exclude a massive body in the asteroidal 
region to be defined as a planet. Therefore, this constraint is automatically included in our first 
constraint. Their second additional constraint is the Earth's water content. They showed that the 
mass supply from asteroidal bodies to terrestrial planets decreases with increasing eccentricity 
of Jupiter. As a result, the water content of an Earth-analog is insufficient in their EJS and 
EEJS simulations. On the other hand, in our EJS and EEJS simulations with the nebular gas, a 
large fraction of mass of an Earth-analog comes from the asteroidal region, while protoplanets 
which formed from planetesimals originally in the terrestrial region spiral into the Sun by Type 
I migration (Sec. 3.1.1). Probably, the problem here is rather that we need to consider how to 
remove excess water from our Earth-analogs. 



3.2.1. Comparison between EJS and CJS simulations 

The properties of the planetary systems from the EJS and CJS simulations are summarized 
in Table 1 . Also, Figures 8 and 9 show the snapshots on the a-e plane at the end of all EJS and [Table 1] 
CJS simulations, respectively. For both EJS and CJS simulations, the total mass and the typical [Figs. 8 
size of planets decrease with increasing gas decay time r decay because of Type I migration. When and 9] 
the surface density is large, protoplanets rapidly grow and then spiral into the Sun. Therefore, 
the final total mass M T final depends on the initial surface density only weakly for a certain r decay . 
This trend is particularly evident when r decay > 2 Myr, because the mass depletion due to Type 
I migration propagates to the outer asteroidal region for large Tdecay The number of planets 7V p 
is usually 3 or 4. In the EJS simulations with r decay = 5 Myr, however, 7V p is smaller; in these 
simulations, the v 5 resonance effectively sweeps most of the remnant mass to its current location 
(0.61AU), and a planet with dominant mass (0.3 - 0.5 M e ) forms there. 

We plot the final planets on the plane of semimajor axis versus mass (Fig. 10). The crosses [Fig. 10] 
on the figure represent the current terrestrial planets. The EJS simulations clearly have better 



2 Smaller fi mp (- 30 Myr) is not ruled out from the terrestrial W/Hf records because the degree of metal-silicate 
equilibration at the Moon-forming impact is still in debate (Jacobsen et al., 2009). In this case, however, one needs 
an alternative interpretation for the prolonged differentiation of the Moon (Touboul et al., 2007), other than the late 
Moon-forming giant impact. 
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radial mass concentration than the CJS simulations. In the CJS simulations, the effect of the 
V5 resonance is negligible so that we inevitably have a massive planet around 2 AU, quite in 
contrast to the EJS simulations. In both EJS and CJS simulations, sizes of outer planets are 
smaller when the initial surface density in the asteroidal region is small, as is the case for small 
Mj or large p. In addition, for the EJS simulations, the radial mass concentration is indirectly 
affected by the initial spatial distribution of planetesimals via their interactions with Jupiter. 
Asteroidal particles reduce ej either by angular momentum exchange at the V5 resonance or by 
energy exchange in close encounters with Jupiter; in both cases, the decrease in gj is given as 
Aej ~ m/(mjej), where nij is Jupiter's mass. With decreasing ej, the effect of the v 5 resonance 
decreases (Nagasawa et al., 2005). If the total mass in the asteroidal region is small due to small 
M T , large p, or large Td eC ay, the decrease in ej is small (see final ej in Table 1). Such systems 
can obtain high radial mass concentration since Jupiter's perturbation remains effective. This 
accounts for the large difference in the radial mass concentration between cases with M T = 5 
and 10 M @ in the EJS simulations (Fig. 10). The v 6 resonance is also important for the radial 
mass concentration, and it also works effectively when e s is large (we will discuss its effect at 
Fig. 17 more in detail). Since Saturn's eccentricity e s is primarily determined by the eccentricity 
forced by Jupiter, e s is usually comparable to or slightly larger than e h 

Figure 1 1 shows the plot of the final planets on the plane of mass versus orbital eccentric- 
ity. The small orbital eccentricities are well reproduced in the EJS simulations with M T = 5 [Fig. 11] 
M ffi but not in those with M T = 10 M 9 . This difference is closely related to the radial mass 
concentration. In the case with M T = 5M®, the V5 resonance effectively sweeps most of the 
bodies and giant impacts between protoplanets complete early while there remains a sufficient 
amount of gas and planetesimals. On the other hand, the systems starting with M T = 10 M @ 
radially expand even after the secular resonance passage and Jupiter continuously perturbs the 
outer planets. These planets eventually trigger orbital instability and cause giant impacts after 
both gas and planetesimals significantly deplete. Giant impacts tend to occur very late in the 
CJS simulations. However, small asteroidal bodies remain for a long time here as well and tend 
to damp eccentricities of the planets. This is why planets in CJS simulations have moderately 
small eccentricities both for M T = 5 and 10 M B . For the CJS simulations, we found a simple 
correlation between the mass weighted mean orbital eccentricity, e m , and the mass weighted 
mean semimajor axis, a m (Fig. 12). The longer the gas decay time, the larger a m , because mass [Fig. 12] 
depletion due to Type I migration propagates outward with time. Since Jupiter perturbs outer 
planets more strongly than inner planets, we have this simple correlation. The CJS simulations 
always produce too large a value for a m , while the EJS simulations well reproduce the current 
value of a m except for runs with r decay = 5 Myr where a m becomes too small. 

Figure 13 shows angular momenta and impact velocities of all potentially Moon-forming 
impacts as functions of time. The most important result is that giant impacts occur most com- [Fig. 13] 
monly at ~ 10 Myr in the EJS simulations whereas at ~ 100 Myr in the CJS simulations. In 
the EJS simulations, some early giant impacts occur between protoplanets trapped in the v 5 
resonance. It is more common, however, that giant impacts occur slightly later when orbital 
eccentricities of planets have sufficiently damped after passage of the resonance so that mutual 
gravitational attraction is enhanced. In the CJS simulations, giant impacts usually occur after 
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gas completely depletes and planetesimals in the asteroidal region also significantly (but not 
completely) deplete. The time of the last Moon-forming impact in each simulation is shown 
in Table 1. In fact, last giant impacts in some of the EJS simulations occur very late. These 
systems, however, inevitably have either large eccentricities or low radial mass concentrations. 

When we see a single simulation in Fig. 13, the impact velocity, i) imp , is as small as the 
mutual escape velocity u esc at early time, but Vi mp /v esc increases with time so that its value for 
the last impact is notably larger than unity. Some simulations, however, have final impacts 
with Vi mp /v esc - 1, which is favorable for the small iron content of the Moon (Canup, 2004, 
2008). When the orbital separation of two planets, normalized by the mutual Hill radius, slowly 
decreases due to their growth or weak perturbation from other planets, an impact between them 
with Vi mp /v esc - 1 can occur after the Jacobi integral in Hill's equation slightly exceeds the 
potential at the Hill's sphere. Orbits of other protoplanets are affected only weakly during the 
time that two colliding protoplanets interact. On the other hand, if three or larger number of 
planets interact strongly, their relative velocity increases to the escape velocity. In this case, the 
typical impact velocity is Vi mp /v esc - V2. Further higher Vi mp /v ssc are possible between small 
protoplanets stirred by larger planets or if Jupiter's perturbation is strong. The former type 
of slow giant impacts usually occur in early stages while the latter type of high- velocity giant 
impacts occur after depletion of both gas and planetesimals (Fig. 13). Systems with S d and e m 
as small as those for the current system do not usually experience the latter type of impacts. 
Therefore, not only the small orbital eccentricities of the Earth and Venus but also the small 
iron content of the Moon are likely to support the presence of some damping forces at the time 
of the Moon-forming impact. 



3.2.2. Cases with larger eccentricities of giant planets: EEJS and CJSECC simulations 

As shown in Table 1, ej (and e s ) decreases via interactions with particles in the asteroidal 
region. Remnant planetesimals orbiting outside of Jupiter's orbit may also reduce e } , although 
we do not include these bodies in our simulations. Therefore, the initial ej and e s are likely to 
have been larger than those adopted in the EJS simulations, if their initial semimajor axes were 
(nearly) the same as the present values. On the other hand, if their initial semimajor axes were 
those given in the Nice model, ej and e s might have been nearly zero throughout the terrestrial 
planet formation, but larger initial ej and e s are not excluded. The summary of simulations with 
larger initial ej and e s , namely the EEJS and CJSECC simulations, is given in Table 2. The [Table 2] 
snapshots of planets in the end of simulations in the a-e plane are also shown in Figs. 14 (the 
EEJS runs) and 15 (the CJSECC runs). [Figs. 14 

and 15] 

In the EEJS simulations, the v 5 resonance is so strong that a m usually becomes too small. 

The number of planets is typically one or two, because protoplanets experience large orbital 
excursions. None of the simulations have S c as large as the present value, because the semimajor 
axis of the innermost planet is usually too small. The innermost planet is pushed toward the Sun 
by outer planets which are swept by the v 5 resonance. The potentially Moon-forming impacts 
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usually occur too early and impact velocities are too high (Fig. 16). Clearly, the initial ej of our [Fig. 16] 
EEJS simulations is too high. If ej is somewhat smaller but still larger than the present value, 
the V5 resonance may work much better (Thommes et al., 2008). 

The initial orbital eccentricities of Jupiter and Saturn in the CJSECC simulations are sim- 
ilar to those in the EJS simulations. The orbital separation between Jupiter and Saturn in the 
CJSECC simulations is narrower than that in the EJS simulations. This makes the v 5 and v 6 
resonances in the CJSECC simulations more distant from the Sun than those in the EJS sim- 
ulations. Because of this, the averaged a m and t[ mp in the CJSECC simulations become larger 
than those in the EJS simulations (Tables 1 and 2). Figure 17 shows a comparison between a 
CJSECC and an EJS simulation with the same parameters except for the orbits of Jupiter and 
Saturn. The v 5 resonance sweeps planetesimals effectively at 3 Myr but no longer works at 5 [Fig. 17] 
Myr in both cases. In the EJS simulation, however, the v^ resonance (~ 2 AU) still sweeps 
planetesimals toward the terrestrial region. On the other hand, the location of the v^ resonance 
in the CJSECC simulation is too far (~ 3 AU) to contribute to radial mass concentration. In 
the CJSECC simulation, accretion of outer bodies, which survived the sweeping by the V5 res- 
onance, proceeds slowly and this makes the final giant impact late. Surprisingly, 7 out of 10 
CJSECC simulations in which potentially Moon-forming impacts occur satisfy the time of the 
last giant impact suggested from isotope records (Table 2). However, iWp/fesc's of most of these 
late impacts are too high (Fig. 16) and also S c 's of these systems are too small (Table 2). 

As we have seen, simulation outcomes strongly depend on Jupiter's eccentricity, ej. Through 
interactions with small bodies, ej decreases. Clearly, the effect of Jupiter's perturbation is more 
important if e } in the late stage is larger regardless of initial values of e } . In Fig. 18, we plot 
the angular momentum deficit S j, the radial mass concentration parameter 5 c /a^, and the time 
of the last one of potentially Moon-forming impacts t{ mp versus the value of ej at the end of a 
simulation, ej,fi n ai (averaged over a few Myr). We find that S c /a 2 n more appropriately represents [Fig. 18] 
the radial mass confinement than does S c because systems consisting of small planets with large 
a tend to obtain large S c while b min is large (particularly in the CJS simulations). Simple linear 
X 2 fits suggest that S c la 2 m and t[ mp increases/decreases with ej^nai, although the scatter of the 
simulation data are very large. The problem is that S c /a 2 n can be as large as the present value 
when ej,finai is large whereas t- imp can be as large as that suggested from isotope records only 
when ej,fi na i is small enough. For Sd, we fit the data with a quartic function as two competing 
effects are expected: large e } means stronger perturbation while quick accretion due to large 
e } means stronger damping forces due to remnant gas and planetesimals. We find that Sd is 
nearly flat for small ej ; fi na i and increases with increasing ej,fi na i when ejfi na i is larger than the 
present value. In an ensemble average sense, S d is usually much larger than the current value, 
although some simulations have S d consistent or even smaller than the present day value. The 
dependence of all three properties on ejfi na i is stronger for the CJS+CJSECC simulations than 
for the EJS+EEJS simulations. It is unclear to us, however, whether this is really physically 
meaningful because fits for the CJS+CJSECC simulations are not very reliable with only a few 
points at large e } . 
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4. Discussion 

None of our simulations succeeded in satisfying all of the three constraints listed in Sec. 3.2. 
The most serious issue we found is a trade-off between the radial mass concentration parameter, 
S c/o-m, and the time of the potentially Moon-forming impact, t- imp (the upper panel of Fig. 19). 
The large S da^ in the current system is due to the small masses of Mars and Mercury and the [Fig. 19] 
small orbital separation between the Earth and Venus. It is hard to satisfy both of these condi- 
tions in a system with large t{ mp . Since the evolution time scale is longer in the outer regions, 
giant impacts can occur late when a large amount of mass stays in the outer region for a long 
time. Such a system usually has outer planets more massive than Mars. Also, if a planetary 
system slowly evolves with weak perturbations from Jupiter (e.g., the CJS simulations), the 
orbital separation between neighboring planets normalized by the mutual Hill radius inevitably 
becomes larger than that between the Earth and Venus (26.3) owing to the mutual orbital re- 
pulsion. If ej is as large as the present value or larger, S c la 2 m can be as large as today's value 
because Jupiter tends to radially compress a system. However, accretion completes too rapidly 
in such a case. 

There is also a trade-off between the angular momentum deficit S 4 and t[ mp (the lower 
panel of Fig. 19). Namely, in order to obtain small Sj, a final large impact needs to occur while 
some damping forces work effectively. Otherwise, 5 ^ ends up with much larger than the current 
value. This trend is particularly clear for the EJS and CJS simulations (filled and open circles). 
Small S d even without any damping forces seems to be possible but such a case is statistically 
rare. 

Our failure indicates that there may be some missing physics in our simulations or flaws 
with the assumptions that we have adopted. In the following, we suggest three possible scenar- 
ios/solutions for future studies. 

In the case with an eccentric Jupiter, a possible solution is that Tdecay is much longer and 
Type I migration is much slower than we assumed. This situation actually corresponds to that 
adopted in the original dynamical shake up model of Nagasawa et al. (2005) and Thommes et al. 
(2008). In their simulations (Td eC ay = 3 and 5 Myr), the giant impacts usually occur at t =* 4 - 5 
Td eca y when the v 5 resonance reaches near 1 AU. Therefore, Td eC ay > 10 Myr is required for the 
late Moon-forming impact indicated from isotope records (t imp > 50 Myr; Touboul et al. 2007). 
If Tdecay = 10 Myr, the gas density is still ~ 1% of the minimum mass solar nebular at 50 Myr, 
and observations suggest that such a gas disk is exceptionally long-living (Pascucci et al. 2006). 
Slightly smaller t[ mp (- 30 Myr) is not excluded from the terrestrial W/Hf records only (Jacobsen 
et al., 2009), but this still requires a long- surviving disk (Tdecay ~ 6 Myr). Recent theories 
indicate significant changes in Type I migration rate due to the non-isothermal effect (Baruteau 
and Masset, 2008; Paardekooper and Mellema, 2008) and the non-linear effect of the co-rotation 
resonance (Masset et al., 2006; Paardekooper and Paparoizou, 2009). These mechanisms cause 
slow down of inward migration or even outward migration if the radial entropy gradient is 
negative for the former mechanism or if the radial surface density gradient is nearly flat or 
positive for the latter one. Therefore, if a local density maximum exists, planetary migration 
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slows down or protoplanets even tend to migrate toward this point (Kretke et al., 2009). A 
local density maximum can exist at the boundary between the optically thin inner region and 
the optically thick outer region (i.e., the inner boundary of the so-called dead-zone; Kretke and 
Lin, 2007; Brauer et al., 2008). Protoplanets might have survived in the protosolar nebular for 
a long time with such mechanisms. 

An alternative possibility for the case with an eccentric Jupiter is that the nebular gas 
inside of Jupiter might have dissipated very rapidly after the formation of Jupiter. This situation 
corresponds to previous /V-body simulations for the late stages without gas (Chambers, 2001; 
O'Brien et al., 2006; Raymond et al., 2009). Without gas, t{ mp tends to be larger than those in 
our simulations, because the effect of V5 is less signifiant. Non-uniform gas dissipation models, 
such as inside-out and gap-opening models discussed in Nagasawa et al. (2000), probably have 
similar outcomes, because v 5 passes a certain location after the gas dissipates there. Raymond 
et al. (2009) showed that some of their EEJS simulations nearly satisfy all the three constraints 
listed in Sec. 3.2. In fact, these simulations satisfy the small mass of Mars, but not the small 
orbital separation between the Earth and Venus. This is why the values of S c in all of their 
simulations are smaller than the current value. This problem might be resolved, however, if 
they adopt even higher resolution with smaller planetesimals (e.g., Morishima et al., 2008) 
and inelastic rebounds for high-speed or grazing collisions (Agnor and Asphau, 2004; Asphau 
et al., 2006). Because these effects are likely to enhance the damping effects and reduce the 
orbital excursions of planets, the orbital separation between neighboring planets will probably 
be smaller. Very quick or non-uniform dissipation of gas seems to be favorable, at least, in the 
asteroidal region in the case with large e } . As we have shown, the v 5 resonance is so strong 
in the asteroidal region, except the CJS simulations, that the radial compositional structure in 
the asteroidal region is destroyed even for r decay = 1 Myr. The gas dissipation of the order of 
0.1 Myr or shorter seems to be necessary for the uniform dissipation model (Nagasawa et al., 
2000). Jupiter might have opened a wide gap around its orbit, extending to the entire asteroidal 
region immediately after its formation (however, Td eca y = 1 Myr or larger may be still possible 
if ej is small, when v 5 passes the asteroid region, and increases later). 

In the case of a circularly orbiting Jupiter, the initial planetesimals needed to be radially 
highly concentrated in order to satisfy a large S c . In other words, the depletion of mass in the 
asteroidal region needed to occur before the formation of planetesimals. Such a concentration of 
planetesimals may be possible if planetesimal formation occurs in radially limited regions after 
the increase of the surface density due to dust migration (Youdin and Shu, 2002; Youdin and 
Chiang, 2004; Kretke and Lin, 2007; Brauer et al, 2008). Morishima et al. (2008) performed 
/V-body simulations starting from highly concentrated disks in the absence of gas, and found 
that some of the resulting planetary systems were similar to the current solar system. The 
problem of their simulations is that accretion usually completes very quickly ( < 10 Myr), as 
well as ourEJS simulations. Nevertheless, one simulation (run 1 of Morishima et al., 2008) has 
a very late giant impact. This system unfortunately has slightly too large an orbital separation 
between the two largest planets. Hansen (2009) also made simulations similar to Morishima et 
al. (2008) and found the mean time of last giant impacts to be 45 Myr. However, the lower limit 
to the impactor mass (0.02 M @ ) used for his definition of giant impacts is probably too small 
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for formation the Moon. Further explorations for planetary accretion starting from radially 
concentrated disks are necessary. 



5. Conclusions 

We have performed iV-body simulations which follow the growth and accretion histories 
of the terrestrial planets, including the effects of the giant planets and the interaction with gas 
in laminar disks. We vary initial total mass and spatial distribution of planetesimals, the time 
scale of dissipation of a gas disk, and the orbits of Jupiter and Saturn. It is a remarkable success 
of numerical simulations, that planetary systems with global properties similar to the terrestrial 
planets can be obtained starting with a proto-planetary disk. We are now at the stage where fine 
details of the models can be compared with solar system constraints. In order to retain a final 
total mass as massive as the present total mass (~ 2 Earth mass) in the presence of standard 
Type I migration, the time scale of dissipation of gas needs to be 1-2 Myr. 

Using the results from a large number of simulations we compared the evolutionary his- 
tories of the resulting planetary systems with the following three properties of the terrestrial 
planets: (1) the high radial mass concentration, (2) the small orbital eccentricities and (3) the 
very late Moon-forming impact (~ 100 Myr). The orbital eccentricities of the planets can be as 
small as their present values if giant impacts finish whilst there is still remnant gas or planetes- 
imals, as expected. The most serious issue we found was a trade off between the radial mass 
concentration and the late Moon-forming impact. If the orbital eccentricities of Jupiter are as 
large as the present value, most of bodies in the asteroidal region are swept up to the terrestrial 
region owing to inward migration of the secular resonance. Although this mechanism helps 
planetary systems to have high radial concentrations in the end, giant impacts between proto- 
planets most commonly occur too early (~10 Myr). On the other hand, if the orbital eccentricity 
of Jupiter is close to zero as suggested in the Nice model, the effect of the secular resonance is 
negligible and a large amount of mass stays for a long period of time in the asteroidal region. 
Although giant impacts usually occur around 100 Myr, we inevitably have an Earth size planet 
at around 2 AU in this case. It is very difficult to obtain spatially concentrated terrestrial planets 
with very late giant impacts, as long as we include all the effects of gas and assume initial disks 
similar to the minimum mass solar nebular. 
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- 


0.004 




5 


2 


10 


3 


0.55 


16.7 


51.7 


1.49 


0.069 


0.0167 


- 


0.003 




5 


2 


5 


4 


0.53 


32.0 


61.7 


1.50 


0.117 


0.0159 


- 


0.003 


ss 








4 


1.98 


89.9 


26.3 


0.90 


0.038 


0.0018 


50-150 


0.048 



Table 1: Simulation parameters and outcomes for the EJS and CJS simulations. Table columns 
are the orbits of Jupiter and Saturn, the gas decay time scale Tdecay, the power-law index for 
the radial gradient of the initial surface density of planetesimals p, the initial total mass of 
planetesimals M T , the number of planets N p , the final total mass of planets M T ,finai> the radial 
mass concentration statistics S c , the mass weighted mean semimajor axis a m , the minimum 
orbital separation between neighboring planets b m [ D normalized by the mutual Hill radius, the 
mass weighted mean eccentricity e m , the angular momentum deficit S^, the time of the last 
potentially Moon-forming impacts t[ mp , and the final orbital eccentricity of Jupiter ej ( fl na i. Any 
orbital elements averaged over a few Myr are used. Each value is boldly marked when 3 < 
N p < 5, 1.0 < M T ,finai/M ffi < 3.0, S c > 45.0, 0.75 AU < a m < 1.05 AU, b min < 35, e m < 0.076, 
Sd < 0.0036, or 50 Myr < P imp < 150 Myr. SS represents the solar system. 



32 



JS orbits 


^"decay 

(Myr) 


P 


(M E ) 


N v 


M T ,final 

(M E ) 


s c 


^min 


(AU) 


C m 


s d 


'imp 

(Myr) 


<?J,final 


EEJS 


1 


1 


10 


3 


3.70 


16.6 


38.1 


0.72 


0.149 


0.0219 


180.0 


0.054 




1 


1 


5 


4 


3.41 


22.6 


30.8 


0.85 


0.053 


0.0023 


62.8 


0.054 




1 


2 


10 


2 


2.66 


40.7 


59.8 


0.86 


0.093 


0.0114 


17.6 


0.060 




1 


2 


5 


1 


1.50 


- 


- 


0.64 


0.089 


0.0113 


26.1 


0.074 




2 


1 


10 


2 


1.71 


36.7 


64.5 


0.56 


0.111 


0.0065 


263.9 


0.067 




2 


1 


5 


4 


2.15 


31.0 


24.8 


0.75 


0.077 


0.0069 


28.3 


0.077 




2 


2 


10 


2 


0.89 


34.0 


154.1 


0.65 


0.289 


0.0427 


23.3 


0.077 




2 


2 


5 


2 


1.23 


52.6 


72.2 


0.95 


0.196 


0.0572 


12.7 


0.046 




3 


1 


10 


2 


0.63 


38.4 


85.1 


0.66 


0.170 


0.0211 


0.6 


0.071 




3 


1 


5 


1 


0.59 


- 


- 


0.72 


0.597 


0.2090 


21.5 


0.074 




3 


2 


10 


2 


0.67 


46.6 


140.0 


0.86 


0.478 


0.1316 


22.6 


0.082 




3 


2 


5 


1 


0.69 


- 


- 


0.73 


0.411 


0.0916 


20.9 


0.089 




5 


1 


10 


1 


0.33 


- 


- 


0.61 


0.171 


0.0159 


- 


0.084 




5 


1 


5 


2 


0.67 


31.1 


49.6 


0.63 


0.061 


0.0144 


- 


0.074 




5 


2 


10 


1 


0.15 


- 


- 


0.72 


0.234 


0.0356 


- 


0.090 




5 


2 


5 


1 


0.21 


- 


- 


0.65 


0.263 


0.0377 


- 


0.093 


CJSECC 


1 


1 


10 


4 


3.72 


27.0 


35.3 


1.28 


0.041 


0.0019 


33.4 


0.004 




1 


1 


5 


3 


2.74 


32.8 


31.7 


1.11 


0.180 


0.0192 


111.9 


0.007 




1 


2 


10 


3 


2.00 


42.6 


48.6 


1.18 


0.077 


0.0076 


50.5 


0.011 




1 


2 


5 


3 


1.96 


50.8 


56.5 


0.95 


0.045 


0.0033 


71.0 


0.014 




2 


1 


10 


5 


2.27 


10.9 


25.0 


0.87 


0.057 


0.0024 


106.3 


0.014 




2 


1 


5 


4 


2.17 


20.7 


33.0 


1.06 


0.065 


0.0038 


110.1 


0.005 




2 


2 


10 


4 


1.43 


49.9 


32.2 


1.13 


0.053 


0.0022 


- 


0.015 




2 


2 


5 


3 


1.46 


59.6 


25.0 


0.89 


0.089 


0.0060 


16.5 


0.014 




3 


1 


10 


3 


1.11 


64.8 


34.5 


1.33 


0.064 


0.0219 


180.4 


0.006 




3 


1 


5 


3 


1.44 


27.3 


65.4 


1.12 


0.076 


0.0060 


80.8 


0.013 




3 


2 


10 


2 


0.71 


89.4 


54.3 


1.02 


0.139 


0.0112 


- 


0.032 




3 


2 


5 


2 


0.95 


85.6 


105.0 


0.93 


0.160 


0.0132 


49.2 


0.032 




5 


1 


10 


3 


0.62 


12.4 


121.7 


0.78 


0.202 


0.0476 


- 


0.013 




5 


1 


5 


3 


0.92 


61.8 


31.6 


0.87 


0.025 


0.0007 


- 


0.037 




5 


2 


10 


2 


0.41 


67.5 


74.1 


1.08 


0.121 


0.0190 


- 


0.027 




5 


2 


5 


2 


0.50 


110.4 


57.9 


1.02 


0.061 


0.0038 


- 


0.029 


ss 








4 


1.98 


89.9 


26.3 


0.90 


0.038 


0.0018 


50-150 


0.048 



Table 2: The same as Table 1, but for the EEJS and CJSECC simulations 



-33- 



io 8 
io 7 

't 10 8 
cC 

o 

^ 10 5 

M 

I 

1000 

100 

10 5 



10 4 



71 

S-. 

o 



'1000 



100 



10 



Migration time 





Damping time 

0.01 
^^" 0.05^ 

^^S ^'0.1 /- 
^-^"0.3 

i i nun i i i innl i i i mill 



Ceres 



• • • 

Moon Mars Earth 



■f-^ 



\ 



V 



111 I 



1L_ 



1Q 1B 1Q19 1Q Z0 1Q 21 10^2 1QZ3 1 Q24 1 Q25 1 QZ6 1 Q27 1 QZB 1 Q29 

Mass m (g) 



Fig. 1 . — Time scales for radial migration T m j g and damping of eccentricity Td am p due to gas drag 
and tidal interaction (r m i g = a[(Ja/J0dra g + (^/^0tid] _1 and Tdamp = e[(de/dt)di ag + (de/dt) t id\~ l ) 
at 1 AU with S gas = 2000 g cm" 2 and i = 0. The migration and damping rates due to gas drag, 
(da/dt) dmg and (de/dt) dmg , are derived after Adachi et al. (1976) with some corrections given in 
Inaba et al. (2001). The labels to the lines represent e. For e = 0.05, we also plot the time scales 
with the enhancement of the gas drag force for m = 0.0025 M@, m\ = 0.01 M 9 , and m G i = 10 19 
g (see Sec. 2.2.1). The masses of solar system bodies are also plotted in the upper panel. 
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Fig. 2. — Snapshots of an EJS simulation (r decay = 1 Myr, p = 2, and M T = 5 M e ) on the a - e 
plane. The size of each body is proportional to the radius, except Jupiter's size is modified to 
the Earth size on this plot. The color represents spin rate to: to > o> cr i t , a» cr it > co > 0.5 o> cr j t , 0.5 
cd cr it > co > 0, and to = for red, blue, green, and black bodies, respectively, where to CT i t is the 
spin rate with which the centrifugal force balances with the gravity at the surface. The long- 
dashed and short-dashed lines are locations of the v 5 and V(, resonances, respectively. The black 
dashed curve is given by a{\ + e) = ay, particles in the right side of this curve experience orbital 
crossing with Jupiter. Note the difference in scales of e between the top panels and others. 
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Fig. 3. — Time evolutions of the total masses of small, big, and all particles (the upper panel) 
and the mass weighted mean eccentricities of small and big particles (the lower panel) for an 
EJS simulation. We define mass of a small (big) particle to be below (above) 2 x 10 26 g. 
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Fig. 4. — The time evolutions of the masses (the upper panel) and the orbital excursions (the 
lower panel) of surviving planets in an EJS simulation. The numbers of planets are given in the 
order of their masses in this figure. 
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Fig. 5. — The same as Fig. 2, but for a CJS simulation. Except for orbits of Jupiter and Saturn, 
other parameters are the same as those in Fig. 2. 
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Fig. 6. — The same as Fig. 3, but for a CJS simulation. 
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Fig. 7. — The same as Fig. 4, but for a CJS simulation. 
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Fig. 8. — Snapshots of all EJS simulations on the a - e plane at the end of simulations. In each 
panel the numbers after t, p, and m represent r decay in units of Myr, p itself, and M T /M e . The 
colors of bodies represent the spin rates as well as Fig. 2. In the panel at the upper left corner, 
we plot the current terrestrial planets with dashed circles. 
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Fig. 9. — The same as Fig. 8, but for the CJS simulations. 
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Fig. 10. — Mass vs. semimajor axis for the EJS and CJS simulations. The filled and open 
symbol are for M T /M m = 5 and 10. The color represents the gas decay time: T decay = 1, 2, 3, and 
5 Myr for black, red, blue, and green symbols. Circles and squares are for p = 1 and 2. Cross 
marks represent the current terrestrial planets. 
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Fig. 11. — Eccentricity vs. mass for the EJS and CJS simulations. See Fig. 10 for symbols. 
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Fig. 12. — Mass weighted mean eccentricity e m vs. mass weighted mean semimajor axis a m for 
the EJS and CJS simulations. The dotted vertical and horizontal lines indicate the values for the 
current terrestrial planets. See Fig. 10 for symbols. 
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Fig. 13. — Angular momenta and velocities of potentially Moon-forming impacts as functions 
of time for the EJS and CJS simulations. See Fig. 10 for symbols. 
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Fig. 14. — The same as Fig. 8, but for the EEJS simulations. 
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Fig. 15. — The same as Fig. 8, but for the CJSECC simulations. 
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Fig. 16. — The same as Fig. 13, but for the EEJS and CJSECC simulations. Some impacts in 
the EEJS simulations have too large L imp (> 6 L EM ) and only their iv p 's are shown. 
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Fig. 17. — Comparison between the EJS and CJSECC simulations. The same parameters are 
used for the two simulations except for orbits of Jupiter and Saturn: Td eC ay = 1 Myr, p = 1, and 
Mj = 5 M m . The format for colors and lines is after Fig. 2. 
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Fig. 18. — Angular momentum deficit Sd (upper panel), radial mass concentration parameter 
S c /af n (middle panel), and time of the last one of potentially Moon-forming impacts ?; mp (lower 
panel) vs. final eccentricity of Jupiter ej^nai- Filled circles, filled squares, open circles, and open 
squares are results from the EJS, EEJS, CJS, and CJSECC simulations, respectively. The color 
represents the gas decay time: Td eC ay = 1>2, and 3 Myr for black, red, and blue symbols. A 
fit to a quartic function is applied to logSj, while linear fits are applied to 5 c /a^, and log? imp 
for sets of the EJS+EEJS and CJS+CJSECC simulations, respectively. Horizontal and vertical 
dotted lines indicate the current values (for t imp , a possible range indicated from isotope records 
is between two dotted lines) 
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Fig. 19. — Radial mass concentration parameter S c /af n (upper panel) and angular momentum 
deficit Sd (lower panel) vs. time of the Moon-forming impact t imp . The data are the same as 
those in Fig. 18. 



